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ABSTRACT 

In  the  context  of  off-road  vehicle  simulations,  deformable  terrain  models  mostly  fall  into  three  categories:  simple  visualization  of  an 
assumed  terrain  deformation,  use  of  empirical  relationships  for  the  deformation,  or  finite/discrete  element  approaches  for  the  terrain.  A 
real-time  vehicle  dynamics  simulation  with  a  physics -based  tire  model  (brush,  beam-based  or  Finite  Element  models)  requires  a 
terrain  model  that  accurately  reflects  the  deformation  and  response  of  the  soil  to  all  possible  inputs  of  the  tire  in  order  to  correctly 
simulate  the  response  of  the  vehicle.  The  real-time  requirement  makes  complex  finite/discrete  element  approaches  unfeasible,  and  the 
use  of  a  ring  or  beam  -based  tire  model  excludes  purely  empirical  terrain  models. 

We  present  the  development  of  a  three-dimensional  vehicle/terrain  interaction  model  which  is  comprised  of  a  tire  and  deformable 
terrain  model  to  be  used  with  a  real-time  vehicle  dynamics  simulator.  The  governing  equations  of  both  models  are  physics-based, 
rather  than  utilizing  popular  terramechanics  models  that  are  empirical.  The  tire  draws  on  a  lumped-mass  model  based  on  a  radial 
spring-damper-mass  distribution.  The  terrain  model  utilizes  Boussinesq  and  Cerruti  soil  mechanics  equations  to  determine  the  pressure 
distribution  and  deformation  of  a  volume  of  soil  as  a  function  of  normal  and  tangent  forces  applied  at  the  soil  surface  by  the  tire.  The 
soil  volume  that  describes  the  terrain  is  discretized  as  a  set  of  vertical  columns  of  soil,  and  the  deformation  of  each  is  modeled  using 
visco-elasto-plastic  compressibility  relationships  that  relate  subsoil  pressures  to  a  change  in  bulk  density  of  the  soil,  which  in  turn 
produces  soil  displacements.  Different  loading  combinations  applied  by  a  tire  passing  over  a  column  of  soil  will  be  reflected  in  the 
state  of  each  volume  of  soil  contained  in  the  column,  rather  than  treating  the  column  of  soil  as  homogeneous  in  the  vertical  direction 
and  only  associating  one  set  of  parameters  with  the  entire  column,  e.g.  a  Bekker  type  model.  Furthermore,  the  time -dependent  elastic 
and  plastic  response  of  the  soil  to  repetitive  compression/rebound  tire  loads  is  also  taken  into  account.  Horizontal  soil 
force/displacement  produced  by  tractive  and  turning  forces  will  also  be  incorporated  into  the  model.  Both  the  vertical  and  horizontal 
force/displacement  relationships  allow  the  calculation  of  total  energy  and  power  required  to  deform  the  terrain.  These  physics-based 
models  will  be  integrated  into  a  real-time  vehicle  dynamics  simulator  and  is  anticipated  to  lead  to  a  realistic  vehicle  dynamic  response 
when  driving  on  off-road,  deformable  terrain  conditions,  especially  when  repeated  loading  occurs  or  non -homogeneous  soil  conditions 
are  present.  Additionally,  the  changes  in  soil  states  can  be  used  to  directly  compute  the  energy  and  power  required  to  deform  the 
terrain.  In  order  to  retain  the  ability  to  run  real-time  simulations,  a  GPU-accelerated  approach  is  considered  to  leverage  the  inherently 
parallel  nature  of  performing  multiple  independent  terrain  geometry  queries  and  soil -mechanics  calculations.  Numerical  experiments 
include  a  single  soil  volume  node  under  a  known  load  and  a  simplified  tire  model  applying  normal  forces  on  the  surface  of  the  terrain. 
Results  are  given  for  the  vertical  plastic  soil  deformation,  and  for  the  power  and  energy  required  to  perform  the  deformations. 

INTRODUCTION 

We  present  the  development  of  vehicle/terrain  interaction  model  to  be  used  with  an  existing  real-time  driving  simulator  at  the  US 
Army's  Tank  Automotive  Research,  Development  and  Engineering  Center  (TARDEC).  The  model  will  make  use  of  existing 
simulation  tools  and  terrain  profile  measurement  information  as  additional  requirements  to  the  development.  The  US  Army  TARDEC 
ride-simulator  (Richmond,  2004)  is  comprised  of  a  multibody  dynamics  vehicle  model  that  is  interfaced  with  an  empirically -based 
deformable  terrain  model.  The  terrain  model  is  able  to  query  lateral  and  longitudinal  location  points  on  a  terrain  database  which 
returns  various  terrain  information  about  a  data  point,  notably  the  undisturbed  terrain  height,  the  change  in  height,  the  Rating  Cone- 
Index  (RCI),  and  the  type  of  terrain,  e.g.  soil  classification  code.  Geometry  associated  with  the  terrain  profile  is  a  combination  of  low- 
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fidelity  terrain  height  measurements  with  superimposed  NURBS  for  high-frequency  content  (Morrison  2004),  and  an  example  terrain 
profile  is  shown  in  Figure  1 .  This  information  is  used  with  empirical  relations  developed  at  the  Engineer,  Research  and  Development 
Center  (ERDC)  to  estimate  the  terrain  sinkage  and  vehicle  tractive  forces  under  a  given  tire  loading  condition.  The  interaction  forces 
at  the  tire/terrain  interface,  which  are  responsible  for  soil  deformation,  are  then  resolved  as  tire-spindle  forces  and  are  passed  back  to 
the  vehicle  model.  This  approach,  based  on  empirical  terramechanics  relations  makes  several  assumptions  such  as,  for  instance,  the 
premise  of  equally  distributed  pressure  under  the  tire  footprint.  The  tire  is  assumed  to  be  rigid,  and  furthermore,  the  soil  mechanics  are 
based  on  static  whole  wheel  approximated  as  a  function  of  vehicle  speed,  which  leads  to  inaccurate  results  for  a  slow  moving  or 
rapidly  accelerating/decelerating  vehicle. 


Figure  1  Terrain  geometry  as  a  heightmap,  after  an  applied  displacement 

The  principle  behind  the  new  Vehicle/Terrain  Interaction  Model  (VTIM)  is  that  the  governing  equations  of  the  model  are  physics- 
based,  instead  of  the  current  empirical  model.  The  inputs  to  the  model  consist  of  vehicle  parameters  and  operating  conditions.  The 
outputs  include  terrain  deformation,  vehicle  power/energy  requirements  and  tire  forces.  The  new  VTIM  builds  around  three 
components:  (i)  a  surface  loading  mechanism  due  to  a  non-rigid  tire,  (ii)  a  rigorous  vertical  soil  stress/strain  integration,  and  (iii) 
lateral  soil  force/displacement  characterization  due  to  soil  surface  shear  and  bulldozing  forces.  They  rely  on  physics-based  tire  and 
terrain  models  and  will  be  integrated  into  the  existing  vehicle  model  and  terrain  database  system. 

The  soil  mechanics  equations  developed  result  in  smooth,  continuous  functions  over  the  affected  volume  of  soil  for  both  vertical  and 
horizontal  soil  deformation;  however,  the  database  will  discretize  these  values  into  a  three  dimensional  grid  under  the  affected  area  for 
database  storage  purposes  and  to  maintain  the  ability  to  run  real-time  simulations.  The  dimension  and  number  of  elements  in  this  grid 
will  be  tailored  to  retain  the  accuracy  of  the  original,  continuous  functions.  Multiple  vertical  and  horizontal  forces  will  be  applied  on 
the  terrain  surface,  depending  on  the  number  of  tire  nodes  in  contact  with  the  terrain.  The  total  horizontal  and  vertical  pressure  for 
each  soil  node  is  determined  by  superimposing  the  contribution  from  each  tire  force  that  causes  subsoil  stress  on  the  particular  soil 
node.  Vertical  sinkage  and  the  soil  states  for  each  soil  node  are  calculated  and  used  to  update  the  terrain  database  for  the  next  time 
step.  They  are  subsequently  used  in  an  energy/power  analysis  that  gauges  the  mobility  and  the  associated  power  requirements  of  the 
vehicle. 

In  the  proposed  approach  real-time  simulations  will  be  possible  due  to  the  inherent  parallel  nature  of  the  tire  node  terrain  queries, 
which  will  utilize  GPU  hardware  to  accelerate  the  computation  of  the  soil  mechanics  response.  The  single-instruction-multiple-data 
programming  model  promoted  by  GPU  accelerators  will  also  be  leveraged  for  tire  forces  and  soil-mechanics  calculations  due  to  the 
node-based  approach  adopted  in  the  proposed  methodology. 

VEHICLE  TERRAIN  INTERACTION  MODEL  COMPONENTS 

This  section  contains  the  physics-based  models  that  will  replace  the  existing  empirical  components  used  in  the  real-time  off-road 
vehicle  dynamics  simulator  at  US  Army  TARDEC.  No  changes  are  made  to  the  vehicle  model,  which  is  already  a  multibody 
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dynamics-based  model.  The  tire  model  is  presented,  followed  by  the  soil  mechanics  model,  which  is  broken  down  into  its  constitutive 
parts  with  an  in-depth  discussion  of  each  component. 

TIRE  MODEL  DESCRIPTION 

In  order  to  enable  real-time  simulation,  the  tire  and  soil  models  are  decoupled.  The  spring-damper  tire  model  combines  a  series  of 
lumped  masses,  linked  to  each  other  and  the  wheel  rim  using  discrete  spring-damper  elements.  The  model  consists  of  two  components, 
which  were  detailed  for  hard  surfaces  in  Negrut  and  Freeman  (1994).  The  first  component  is  the  model  shown  in  Figure  2,  which 
incorporates  both  radial  springs  and  dampers,  and  interradial  springs.  The  second  component  is  the  tire  tread  discretized  into  a  number 
of  lumped  masses.  Lumped  masses  are  connected  to  each  other  and  to  the  wheel  rim  using  tangential  and  lateral  spring  elements.  The 
numerical  results  presented  in  this  work  do  not  contain  tire  dynamics,  thus  the  lumped  mass  tire  tread  description  is  excluded,  but  can 
be  found  in  the  literature. 


Figure  2  Radial-interradial  spring-damper  tire  model 


Each  of  the  spring-damper  elements  composing  the  rolling  radial-interradial  model  can  either  be  an  active  element,  in  contact  with  the 
ground  and  queried  using  the  terrain  geometry  database,  or  a  passive  element.  In  order  to  determine  the  active  elements,  the  terrain 


geometry  is  queried  for  each  lumped  mass  coordinate  tuple,  r  =  {x .  yt  z. } 
the  length  of  the  spring-damper  element  is  calculated  as 


j(r'—r  )7 

(r'—r  ) 

sj\  i  c) 

V  i  c) 

If  there  is  penetration  of  the  terrain  geometry,  then 


(1) 


where  r\  is  the  coordinate  of  Pj  ,  the  ith  terrain  contact  point,  and  r  is  the  coordinate  of  the  center  of  the  wheel,  both  expressed  in 
global  coordinates.  For  each  active  element,  a  Boolean  activity  array  will  be  set  to  true  (1),  while  for  each  passive  element,  the  activity 
array  will  be  set  to  false  (0).  Based  on  wheel  rotational  speed,  the  array  will  be  used  to  minimize  the  number  of  terrain  queries  at  each 
timestep. 

The  radial  force  acting  on  the  ith  mass  element  can  then  be  expressed  as 

Fi=k{R-$i)  +  q(£i+1  - £. )  +  q(£t_x  ~ £ , ) - ci.  (2) 

where  k  and  q  are  the  stiffness  of  the  radial  and  interradial  springs,  C  is  the  radial  damping,  and  R  is  the  undeformed  radius  of  the 
tire.  For  the  undeformed  elements,  the  force  F{  is  zero. 
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SOIL  MECHANICS  MODELS 


The  physics-based  soil  mechanics  model,  developed  by  Ayers  and  Bozdech,  is  a  combination  of  many  different  models  and  is  the 
most  substantial  theory  contribution  to  the  project.  As  such,  this  section  is  broken  into  its  constitutive  parts,  which  include:  vertical 
stress  and  sinkage,  repeated  loading  effects,  shear  stress/shear  displacement  effects.  Power  and  energy  calculations  for  both  vertical 
and  horizontal  soil  deformations  are  also  included. 

Vertical  Visco-Elasto-Plastic  Model 


As  the  soil  surface  contact  area  and  the  vehicle  normal  surface  loads  are  known,  the  subsurface  stress  distribution  needs  to  be 
determined.  The  method  for  calculating  the  stress  in  a  semi-infinite,  homogeneous,  isotropic,  elastic  medium  subject  to  a  vertical 
point  load  applied  on  the  surface  was  first  developed  by  Boussinesq  (Wong,  2001): 

3  i  W 

+  (r/z)2f!r  <3> 

where  <Jz  is  the  stress  at  the  depth  z  ,  W  is  the  point  load,  and  T  is  the  horizontal  distance  from  the  point  load  to  the  point  in 
equation.  When  the  Boussinesq  model  is  applied  for  soil,  there  is  a  tendency  for  the  compressive  stress  in  the  soil  to  concentrate 
around  the  loading  axis.  This  tendency  becomes  greater  when  the  soil  becomes  more  plastic  due  to  increased  moisture  content  or  when 
the  soil  is  less  cohesive,  such  as  sand.  In  view  of  this  phenomenon,  Frohlich  introduced  a  concentration  factor  V  to  Boussinesq’ s 
equations,  for  hard  soil,  V  =  3 ,  for  normal  soil  V  =  4 ,  and  for  soft  soil,  V  =  5  (Koolen  and  Kuipers,  1983).  Other  literature  employs 
concentration  factors  for  hard  soil,  V  =  4 ,  for  normal  soil,  V  =  5  ,  and  for  soft  soil,  V  =  6  (Ayers,  1991;  Wong,  2001).  The 
concentration  factor  has  been  assumed  to  vary  with  the  current  soil  density,  the  minimum  density,  and  the  maximum  density  of  the 
soil  material  (Binger  and  Wells,  1989).  An  equation  was  developed  to  calculate  the  value  of  concentration  factor. 


v  =  6  +  (6-4)(P-Pmin)/(Pmin-Pmax)  (4) 

where  V  is  the  concentration  factor,  P  is  the  current  dry  bulk  density  of  the  soil  (g/cm3),  Prwsi  is  the  maximum  dry  bulk  density  of 
the  soil  (g/cm3),  and  /Vn  is  the  minimum  dry  bulk  density  of  the  soil  (g/cm3).  Introducing  the  concentration  factor,  the  vertical  stress 
under  the  plate  can  be  expressed  in  Cartesian  coordinates  as  (Ayers  and  Riper,  1991): 

vWzv 

°l~  2  n{r2+z2){vl2+l)  (5) 

After  the  subsurface  stress  distribution  is  determined,  soil  strain  (or  compaction)  is  based  on  soil  visco-elasto-plastic  relationships. 
These  relationships  are  to  be  defined  for  the  USCS  terrain  described  in  the  terrain  database.  The  constitutive  (stress/strain) 
relationships  for  these  soils  are  needed  to  define  the  soil  response  (sinkage)  due  to  loading. 

If  bulk  density,  p,  is  plotted  vs.  logi0  applied  stress,  <Ja  ,  the  relationship  will  be  linear.  During  the  compression  procedure,  soil 

particles  will  be  rearranged  and  brought  closer  together.  The  straight  line  is  also  called  the  virgin  compression  curve,  VCC,  and  the 
slope  of  that  line  is  called  the  compression  index,  C  (Larson  et  al.,  1980).  If  the  soil  has  previously  been  subjected  to  a  stress  and  a 
further  stress  is  applied,  the  bulk  density- stress  relationship  will  increase  along  a  line  with  a  slight  slope  (secondary  compression)  until 

it  joins  with  the  VCC  curve.  The  linear  portion  of  the  VCC  in  Figure  3  can  be  described  by  (Larson  et  al.,  1980) 

p  =  pk  +  C \ogU)  (a Jak),  (6) 


where  p  is  the  computed  bulk  density,  pk  is  the  density  at  a  known  stress  <Jk,  <Ja  is  the  applied  normal  stress,  and  C  is  the  slope  of 
the  line.  If  the  VCC  represents  a  dry  soil,  lines  for  soils  at  water  contents  less  than  saturation  will  be  shifted  to  the  left  and  parallel  to 
the  dry  soil  line.  If  the  relationship  in  Eq.  (6)  is  known  for  a  given  water  content,  or  degree  of  water  saturation,  S^,  curves  for  other 
water  contents  can  be  computed  from  the  following  relationship  (Larson  et  al.,  1980): 
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P  =  [pk+ST(S{- sk )]  +  Clog10  (aa  /  <jk ) 


(7) 


where  ST  is  the  slope  of  the  bulk  density  vs.  degree  of  water  saturation  curve  at  a  given  stress,  and  Sj  is  the  desired  degree  of 
saturation.  By  predicting  the  bulk  density  variation  due  to  repetitive  normal  loading,  the  soil  element's  vertical  displacement  can  be 
determined  along  with  the  related  net  energy  and  power  that  must  be  applied  to  the  mass.  A  soil  element  with  a  given  height  ( H)  has  a 
vertical  displacement  (z)  expressed  by  the  following  equation: 


z=(l-(pj  Pl)\H 


(8) 


where  is  the  initial  bulk  density  and  ^is  the  final  bulk  density. 


Power  and  Energy  Calculations 

Positive  soil  element  vertical  displacement  occurs  during  compression  when  px  exceeds  p0  in  Eq.  (8).  The  required  energy  ( E )  to 
displace  the  soil  element  vertically  ( Az )  is  calculated  by  the  following  equation: 

E  =  F  •  Az  =  (cr  •  A)  •  Az  (9) 

where  a  is  the  sub-surface  normal  stress,  A  is  the  soil  element  cross-sectional  area  in  the  horizontal  plane. 

Equation  (9)  is  only  applied  in  the  model  during  compression  (i.e.  positive  Az).  E  is  assumed  to  remain  constant  during  rebound 
because  the  energy  lost  by  the  soil  element  is  assumed  to  be  negligible  in  the  model.  The  power  required  during  compression  of  the 
soil  element  in  a  given  time  interval  (At)  is  determined  by  the  following: 

P  =  E/At  (10) 

P  is  nonzero  only  when  compression  of  the  soil  element  occurs. 

Repeated  Vertical  Loading 

The  time-dependent  elastic  and  plastic  properties  of  the  soil  are  needed  to  fully  describe  the  soil  response  to  tire/track  loads.  Rebound 
properties  are  defined  by  the  soil  rebound  index  or  measure  of  the  decrease  in  density  as  the  applied  load  is  released.  The  time- 
dependent  or  viscous  soil  properties  can  be  defined  using  a  time  constant  of  soil  deformation. 

The  compression  and  rebound  time  constants,  T c  and  Tr ,  must  be  defined  for  the  soil  element  and  the  values  are  indicators  of  the 
time  required  for  the  dependent  variable  p  to  change  from  the  initial  to  the  final  value  during  compression  and  rebound.  The 
following  is  the  general  form  of  the  time-dependent  compression  and  rebound  factor  of  the  model,  respectively: 

(11) 

During  compression,  the  change  in  bulk  density  is  calculated  by  multiplying  the  theoretical  change  in  bulk  density  determined  from 
Eq.  (7)  by  the  time-dependent  compression  factor  from  the  above  expression.  By  doing  so,  the  effect  of  the  duration  that  a  given 
normal  load  being  exerted  on  the  soil  element  is  accounted  for  in  the  model.  If  a  vehicle  travels  at  a  slower  rate,  the  normal  load  is 
applied  over  a  longer  period  of  time;  thus,  the  bulk  density  increase  should  be  greater  at  a  lower  travel  speeds.  The  time-dependent 
factor  attempts  to  account  for  the  time -dependent  nature  of  a  soil  element's  response  to  an  applied  normal  load. 

The  rebound  constant  must  also  be  predicted  for  the  soil  element.  The  value  of  the  rebound  constant  is  the  amount  of  rebound  that 
occurs  when  the  calculated  bulk  density  associated  with  the  applied  normal  stress  is  less  than  the  current  density.  If  the  bulk  density 
increase  is  less  than  the  rebound  constant  multiplied  by  two,  the  magnitude  of  the  rebound  in  the  model  is  estimated  as  half  of  the  bulk 
density  increase  during  compression.  The  rebound  that  occurs  in  a  given  time  interval  is  also  calculated  by  multiplying  the  time- 
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dependent  factor  from  above  by  the  net  rebound  where  t  is  equal  to  zero  when  compression  of  the  soil  element  ends  and  rebound 
begins. 

Lateral/Longitudinal  Shear  Stress/Deformation  Model 

To  produce  mobility,  tractive  tire  forces  must  be  supported  by  the  soil.  These  forces  produce  wheel  or  track  slip  and  produce 
additional  soil  deformation.  The  ability  of  the  soil  to  support  tractive  forces  can  best  be  described  by  both  the  soil  shear  strength 
parameters  (cohesion  and  friction  angle)  and  the  soil  shear  deformation  modulus  (describes  the  shear  stress/deformation  relationship 
of  the  soil).Both  longitudinal  and  lateral  horizontal  shear  forces  are  known  to  produce  vertical  subsurface  stresses,  resulting  in 
increased  vehicle  sinkage,  and  the  phenomenon  is  known  as  “slip  sinkage”.  Field  and  laboratory  tests  of  turning  vehicles  have  shown 
that  lateral  stresses  produce  lateral  displacement  and  increased  sinkage.  This  additional  sinkage  is  considered  “slide  sinkage”. 

The  method  for  calculating  the  vertical  stress  in  a  semi -infinite,  homogeneous,  isotropic,  elastic  medium  subject  to  a  horizontal  point 
load  applied  on  the  surface  was  first  developed  by  Cerruti  (Feda,  1978): 


3  r(cos  0)  H 
2tc\Y  +  (r/z)2]5/27 


where  Gz  is  the  vertical  stress  at  the  depth  z  ,  H  is  the  horizontal  force,  0  is  the  angle  between  the  applied  force  direction  and  the 
direction  to  the  location  of  the  vertical  stress,  and  r  is  the  horizontal  distance  from  the  point  load  to  the  location  of  the  vertical  stress. 
These  horizontal  load  induced  vertical  stresses  will  be  added  to  the  vertical  stresses  produced  from  the  tire/track  vertical  soil  surface 
forces  to  produce  the  total  vertical  stress  at  a  subsurface  point. 

Turning  vehicles  add  even  more  complexity  to  the  vehicle/terrain  interaction  model.  For  wheeled  vehicles,  the  forces  include  the  static 
weight,  plus  any  dynamic  weight  (due  to  weight  shift),  plus  increased  longitudinal  forces  (vehicle  rolling  resistance  increases  with 
turning)  producing  wheel  slip,  plus  centrifugal  (lateral)  forces  due  to  the  turning  radius  and  velocity,  plus  any  forces  produced  as  non¬ 
turning  wheels  or  tracks  are  dragged  across  the  soil  surface.  The  lateral  movement  of  the  wheels  of  a  turning  vehicle,  and  the  slip 
occurring  in  the  longitudinal  direction  will  produce  shear  stress,  and  cause  shear  displacement  of  the  soil.  A  study  of  the  relationship 
between  shear  stress  and  shear  displacement  of  soils  is  necessary  to  understand  both  the  terrain  deformation  and  the  required  energy 
and  power.  Wong  (2001)  indicated  that  there  are  three  types  of  relationship  between  shear  stress  and  shear  displacement.  The  variation 
of  strain-stress  relationship  was  caused  by  the  soil  type,  moisture  content,  and  bulk  density  of  soils.  A  physics-based  shear  stress- 
displacement  model  was  proposed  by  Janosi  and  Hanamoto  (Wong,  2001)  : 

r  =  rma x-(l-e~i/K)  =  (c  +  cr-tmi(/>)-(l-e-j,K)  (13) 

where  T  is  the  shear  stress,  j  is  the  shear  displacement,  C  is  the  internal  cohesion  of  the  soil,  (ft  is  the  angle  of  internal  friction  of  the 
soil  and  K  is  defined  as  the  shear  deformation  modulus.  For  turning  wheeled  vehicles,  the  lateral  shear  stress  resulting  from  the 
centrifugal  force  produces  an  applied  shear  displacement.  The  shear  force  multiplied  by  the  shear  displacement  results  in  the  shearing 
energy  for  the  turning  vehicle.  For  a  turning  tracked  vehicle,  the  track  slides  back  and  forth  a  significant  amount.  This  track  slide  will 
produce  lateral  soil  displacement.  Equation  (13)  can  then  be  applied  to  determine  the  resulting  lateral  forces.  The  power  can  be 
calculated  from  the  force,  displacement  and  time  in  the  same  way  as  the  vertical  power  was  defined  in  Eq.  (10). 

There  is  a  tremendous  loss  of  energy  to  produce  large  amounts  of  soil  deformation  and  relocation.  In  addition,  during  a  turn  the 
outside  wheel  slips  and  the  inside  wheel  is  being  dragged  forward;  thus  the  power  required  for  operating  a  turning  wheeled  or  tracked 
vehicle  is  significant.  The  authors  have  conducted  extensive  research  involving  the  rutting  (both  depth  and  width)  resulting  from 
turning  wheeled  and  tracked  vehicles.  Different  ruts  are  formed  by  the  inside  and  outside  track/tire.  Soil  piling  during  turning  can 
affect  the  vertical  profile  of  the  terrain.  In  addition  to  the  shear  displacement  resulting  at  the  track  shoe/soil  interface,  any  sinkage  will 
create  a  bulldozing  effect  that  pushes  the  soil  to  form  outside  piles  on  both  sides  of  the  tire/track.  The  bulldozing  forces  can  be 
modeled  using  the  Passive  Lateral  Earth  Pressure  calculations,  as  shown  in  Figure  3.  The  forces  resisting  bulldozing  are  the  weight  of 
the  soil  and  the  shear  forces  along  a  defined  failure  plan.  To  simplify  calculations,  negligible  friction  between  the  soil  and  the 
tire/track  interface  can  be  assumed,  resulting  in  a  straight  failure  plane  between  the  bottom  of  the  tire/track  and  the  soil  surface.  The 
angle  of  failure  plane  is  45  -  (f)  12. 
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Figure  3  Force/deformation  relationship  on  a  turning  wheel  or  track  due  to  bulldozing 
The  total  force  at  failure  of  a  sunken  tire/track  shearing  and  displacing  a  soil  wedge  is  defined  as 

F  =  b(y2yZ2N(p+2cZjNi ~)  (14) 

Where  F  is  the  bulldozing  force  at  failure,  y  is  the  soil  unit  weight,  b  is  the  width  of  the  element,  and 


Nv  =tan2(45  +  (z>/2)  (15) 

So  to  determine  the  bulldozing  force  at  soil  failure,  the  soil  strength  parameters,  soil  density,  sinkage  and  width  of  contact  are  needed. 

-3/ 

If  the  bulldozing  lateral  displacement  is  known,  then  the  bulldozing  force  is  reduced  by  (l  —  e  /k),  where  j  is  the  shear  displacement 
and  k  is  the  shear  deformation  modulus.  It  should  be  noted  that  as  the  tire/track  sinks,  Z  increases  and  the  bulldozing  force  also 
increases. 


IMPLEMENTATION  OF  SOIL  MODELS 

At  each  time  step,  the  state  variables  of  the  vehicle  model  will  be  passed  to  the  VTIM,  where  they  will  be  used  in  the  tire  model  to 
determine  the  points  of  contact  with  the  terrain  surface.  Once  the  kinematic  contact  points  are  determined,  the  tire  model  radial  forces 
can  be  determined  using  Eq.  (2),  and  the  points  in  contact  with  the  terrain  can  be  decomposed  into  normal  and  tangential  forces  by 
utilizing  the  terrain  surface  normal  information  contained  in  the  database.  Each  force  can  be  used  to  find  a  subsoil  pressure 
distribution  by  using  the  modified  Boussinesq  and  Cerruti  models,  Eq.  (5)  and  (12),  respectively.  For  A  tire-terrain  contact  nodes,  a 
total  of  2N  pressure  volumes  will  be  calculated,  one  each  for  the  normal  and  tangential  force  components. 

The  tangential  forces  produced  by  slip  between  the  tire  and  the  terrain  are  calculated  on  the  surface  of  the  soil  using  Eq.  (13).  If  the 
vehicle  is  turning,  an  additional  tangential  force  from  bulldozing  is  calculated  using  Eq.  (14).  Combining  the  calculated  tangential 
forces  with  the  displacement  of  the  soil  allows  for  the  calculation  of  energy/power  required  for  the  horizontal  soil  displacement.  The 
energy  and  power  required  for  both  the  vertical  and  horizontal  soil  deformation  constitute  the  total  extra  energy  required  by  the 
vehicle  to  operate  on  the  deformable  soil.  A  flowchart  of  the  calculations  performed  for  each  time  step  in  the  VTIM  is  shown  in  Figure 
4. 

A  given  point  in  the  soil  volume  that  is  affected  by  the  tire  contact  pressure  will  have  a  total  vertical  pressure  applied  to  it  by 
superimposing  the  contributions  of  all  the  vertical  subsoil  pressures  affecting  that  point,  according  to  Boussinesq  and  Cerruti 
equations.  The  soil  in  the  terrain  volume  is  discretized  into  a  three-dimensional  grid  of  rectangles,  where  the  horizontal  dimensions  are 
a  function  of  the  underlying  geometry  query  routines  and  the  depth  is  a  tunable  parameter  depending  on  the  desired  fidelity  of  the  soil 
mechanics  calculations.  Thus  the  soil  point  becomes  a  volume,  and  the  vertical  pressure  applied  at  the  top  of  the  element  can  be  used 
in  conjunction  with  the  current  state  of  the  volume  of  soil  to  calculate  the  change  in  soil  density,  Eq.  (7),  and  vertical  sinkage,  Eq.  (8). 
If  the  soil  volume  is  currently  undergoing  compression,  then  the  required  energy  and  power  are  calculated  with  Eqs.  (9)  and  (10).  If 
the  soil  volume  is  currently  undergoing  rebound,  only  the  change  in  vertical  sinkage  and  soil  density  are  computed. 

The  number  of  calculations  required  to  determine  and  update  the  soil  node  states  is  large  and  can  be  determined  as  follows.  Assuming 
a  predetermined  volume,  V  =  L*W  *  D ,  with  length  L,  width  W  and  depth  D,  of  terrain  nodes  that  includes  the  soil  that  will  be 
affected  by  the  tire  loading  (i.e.,  an  area  slightly  larger  than  the  tire  footprint,  and  a  depth  of  no  more  than  2  feet) ,  a  total  of  2 NV,  or 
2 NLWD  summations  are  required  to  determine  the  vertical  pressure  on  each  terrain  node  affected  by  the  surface  load.  Furthermore, 
LWD  soil-mechanics  calculations  are  required.  Taking  into  account  that  the  spacing  between  the  terrain  nodes  that  describe  the  surface 
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is  small  (to  take  into  account  high-frequency  content  of  the  terrain  profile),  a  huge  number  of  calculations  are  required  to  update  the 
soil  states  due  to  a  tire  load  at  each  timestep. 
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Figure  4  Calculation  Flowchart 


OPPORTUNITIES  FOR  PARALLELISM 

Due  to  the  nature  of  the  equations  above,  it  can  be  seen  that  the  terrain  geometry  is  spatially  independent  -  that  is,  the  terrain  in  one 
area  does  not  depend  upon  that  in  another.  The  same  is  true  for  the  individual  tire  elements.  This  means  that  calculations  on  different 
areas  of  the  terrain  and  different  elements  of  the  tire  can  be  embarrassingly  parallel  -  it  does  not  matter  when  (or  where)  the 
calculations  for  each  area  are  performed,  as  long  as  they  are  all  finished  before  moving  on  to  the  next  set  of  calculations.  This  type  of 
Single  Instruction,  Multiple  Data  (SIMD)  problem  maps  perfectly  to  the  GPU  (Graphics  Processing  Unit)  which  can  typically  have 
hundreds  of  processors  doing  the  same  calculations  on  different  sets  of  data  in  parallel. 

For  example,  given  a  12  inch  wide,  36  inch  diameter  tire,  and  assuming  a  0.2  inch  nodal  resolution  in  each  direction,  there  can  be  up 
to  10,800  nodes  for  which  the  tire  model  needs  to  be  directly  applied  to.  Additionally,  there  are  at  least  that  number  of  nodes  on  the 
terrain  for  which  the  soil  model  must  be  applied  (depending  on  how  many  layers  of  terrain  nodes  are  used).  At  each  time  step  in  the 
simulation,  none  of  the  calculations  within  each  model  depend  on  each  other,  allowing  them  to  be  performed  in  parallel. 

As  mentioned  above,  SIMD  problems  map  perfectly  to  the  GPU  in  that  the  GPU  is  most  efficient  when  all  of  its  hundreds  of 
processors  are  working  on  the  same  calculation,  albeit  using  different  data.  However,  one  issue  with  the  GPU  (and  most  types  of 
processors,  for  that  matter)  is  that  it  can  quickly  become  data  starved  -  it  processes  the  data  much  faster  than  it  can  send  and  receive  it. 
Nvidia  GPUs  are  designed  in  such  a  way  so  as  to  side  step  this  issue.  They  have  a  dedicated  copy  engine  for  moving  data  to  and  from 
memory,  thereby  allowing  the  main  processors  to  focus  on  calculations  and  not  on  moving  data.  As  long  as  the  main  processors  are 
sufficiently  over-subscribed  (that  is,  there  are  thousands  or  tens  of  thousands  of  calculations/threads  waiting  to  be  performed),  the 
latency  of  copying  data  will  be  hidden  as  the  processors  are  never  idle  while  waiting  for  new  data.  While  this  architecture  can  (and  is) 
be  applied  to  CPUs  as  well,  the  GPU  still  has  the  advantage  in  that  its  memory  bandwidth  is  approximately  five  times  that  of  a  typical 
CPU  (120GB/s  for  a  Nvidia  Tesla  C2070  compared  to  20GB/s  for  an  Intel  Xeon  E5520).  Figure  5  shows  approximate  bandwidth 
speeds  for  moving  data  between  different  parts  of  a  computer.  It  can  also  be  seen  that  it  is  desirable  to  move  as  little  data  as  possible 
between  the  CPU  and  GPU  (co-processor)  due  to  the  inherent  latencies  of  the  PCI-e  2.0  system  bus  which  connects  the  system 
memory  on  the  CPU  to  the  device  memory  on  the  GPU. 
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Figure  5  Bandwidth  speed  comparison 

The  MapReduce  framework/methodology  is  used  to  process  this  large  amount  of  data.  MapReduce  is  a  software  methodology 
introduced  by  Google  and  inspired  by  the  map  and  reduce  functions  from  funtional  programming  that  defines  the  workflow  for 
processing  large  amounts  of  data  quickly  using  multiple  processors  working  in  parallel.  The  map  section  of  the  method  splits  the  input 
data  set  so  that  it  is  evenly  distributed  across  all  the  processors.  Each  processor  then  applies  any  user-defined  functions  to  the  dataset 
in  parallel  and  then  passes  the  results  back  into  memory.  The  reduce  section  of  the  algorithm  is  then  applied  to  this  set  of  results  in 
order  to  get  the  final  desired  result.  This  section  could,  for  example,  add  up  all  values  calculated  previously,  sort  them,  and/or  find  all 
unique  values. 

One  of  the  most  common  functions  performed  by  the  VTI  is  querying  the  data  stored  for  multiple  points  of  the  terrain.  In  order  to 
reduce  the  latency  for  copying  data  from  memory  to  the  processors,  multiple  points  are  grouped  together  and  copied  to  the  processors 
in  one  large  chunk.  This  process  is  performed  as  follows:  first,  the  list  of  all  points  required  is  assembled  on  the  GPU  and  evenly 
distributed  to  each  GPU  processor.  Each  processor  will  then  apply  a  basic  spatial  hashing  function  to  each  element.  A  reduction 
algorithm  is  then  applied  to  these  hashes  to  remove  any  duplicates.  All  the  data  for  each  of  these  unique  hashes  is  then  copied  to  the 
processors  for  use  in  the  next  step.  It  should  be  noted  that  this  process  will  bring  in  data  for  points  that  were  not  directly  queried.  This 
is  actually  desirable  because  there  is  a  strong  likelihood  that  the  extra  data  for  those  points  will  soon  be  required  for  future  time  steps. 
Thanks  to  the  multiple  layers  of  cache  in  the  GPU,  this  data  will  then  be  readily  available  to  the  processors  once  they  do  need  it, 
reducing  or  even  eliminating  the  cost  to  copy  it  in  from  memory. 

Once  the  terrain  data  has  been  loaded  on  to  the  GPU  processors,  the  terramechanics  functions  can  be  applied.  Again,  since  these 
functions  may  be  applied  on  an  independent,  point-by-point  basis,  they  may  be  performed  in  parallel  on  each  of  the  GPU's  scalar 
processors. 


NUMERICAL  EXPERIMENTS 

An  example  of  a  test  load  being  repeatedly  applied  to  a  single  soil  volume  is  presented.  This  example  tests  the  response  of  the  soil 
model  to  a  directly  applied  load,  and  thus  the  calculation  of  subsoil  stresses  with  the  Cerruti  and  Boussinesq  equations  is  not 
necessary.  The  applied  normal  stress  to  the  soil  node  was  calculated  using  data  of  a  Stryker,  operating  in  a  CL  soil  at  Ethan  Allen 
Firing  Range,  traveling  at  a  constant  rate  of  1.5,  3.0,  and  6.0  m/s.  As  can  be  seen  by  Figure  6,  the  slower  vehicle  speeds  increases  the 
bulk  density  of  the  volume  of  soil  at  a  greater  rate  initially;  however,  due  to  repeated  loading,  the  rate  of  bulk  density  increase  is 
lowered  due  to  the  increased  density.  This  change  of  rate  of  increase  is  less  apparent  in  the  plot  with  the  higher  vehicle  speed. 
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Figure  6  Bulk  density  and  Applied  stress  as  a  function  of  time,  1.5  m/s  (left)  and  3.0  m/s  (right)  vehicle  speed 

The  total  vertical  displacement,  energy  and  peak  power  due  to  the  deformation  of  the  load  applied  by  the  Stryker  moving  at  1.5  m/s  is 
shown  in  Figure  7.  Notice  that  due  to  the  effects  of  repeated  loading,  a  majority  of  the  energy  required  to  deform  the  soil  is 
accumulated  from  the  first  two  applied  loads.  The  total  vertical  displacement,  energy  and  peak  power  attained  for  each  of  the  three 
vehicle  speeds  is  shown  in  Table  1.  Although  the  peak  power  is  nearly  the  same  of  all  three  vehicle  speeds,  the  total  energy  and 
vertical  displacement  decrease  as  the  vehicle  speed  increases.  This  is  due  to  the  reduced  time  that  the  vehicle  load  is  applied,  and  also 
a  result  of  the  change  in  the  rate  of  bulk  density  increase. 
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Figure  7  Vertical  displacement,  Power  and  Energy  from  vehicle  speed  of  1.5  m/s 

In  order  to  test  the  effect  that  the  Boussinesq  and  Cerruti  equations  have  when  they  are  applied  to  the  soil  database,  a  tire  normal  force 
is  applied  which  is  the  result  of  an  applied  vertical  displacement  on  the  wheel  CG.  A  simplified  tire  model  consisting  of  a  set  of  radial 
springs  applies  a  normal  force  on  the  surface  of  the  terrain.  Also,  the  following  settings  were  used  for  the  terrain  mechanics  settings  of 
all  soil  nodes:  an  assumed  Frolich  stress  concentration  parameter  of  3,  an  initial  bulk  density  of  0.0451591 1  [lb/in3],  and  ak  =  20.356 
[psi]  was  used  for  all  soil  nodes. 

Table  1  Total  vertical  displacement,  Energy  and  Peak  Power  from  various  vehicle  speeds 


Travel 

Speed 

[m/s] 

Correct  Total 
Element  Vertical 
Displacement  [cm] 

Total 
Energy  [J] 

Peak  Power 
[W] 

1.5 

0.095519 

0.039816 

0.180938 

3.0 

0.060483 

0.031553 

0.184590 
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6.0 

0.033694 

0.019222 

0.186478 

The  following  results  are  from  the  test  simulation  detailed  in  the  above  section,  and  are  presented  in  the  following  order:  soil  column 
total  vertical  deformation,  soil  column  energy  to  perform  the  deformation,  and  terrain  surface  profile.  The  maximum  applied  tire  force 
and  subsoil  stress  are  listed  in  table  1  below,  as  well  as  the  energy  required  to  deform  the  soil  node  with  the  maximum  subsoil  stress 
value. 


Table  2.  Max  tire  force  and  subsoil  stress  at  various  steps 


Step  # 

Max  Tire 
Force  [lb] 

Max  subsoil  node 
Stress  [psi] 

Max  Energy 
[lb-in] 

0 

621.7 

505.8 

6.93E-6 

1 

874.4 

763.5 

1.16E-5 

2 

1179.67 

1134.1 

1.89E-5 

3 

1559.1 

1723.4 

3.13E-5 

4 

1996.8 

2476.9 

4.83E-5 

5 

2476.9 

3378.3 

6.89E-5 

The  following  plots  show  the  sum  of  the  deformation  of  the  soil  columns  that  are  affected  by  the  applied  tire  force.  Figure  10  show  the 
time  evolution  of  the  terrain  due  to  the  applied  tire  displacement.  Figure  8  shows  an  isometric  view  of  the  accumulated  vertical  soil 
deformations.  Notice  that  even  though  very  large  tire  forces  and  subsoil  stresses  are  applied  in  this  test  simulation,  the  time  constant 
effect  on  the  compressive  loads  prevents  the  soil  deformation  from  instantaneously  reaching  the  applied  final  deformation  of  1”, 
which  would  be  seen  in  all  the  locations  in  Figure  1  near  the  contact  patch.  This  is  not  the  case,  as  shown  in  Figure  8  where  the 
maximum  soil  deformation  is  less  than  0.9”. 


o 

-0.1 
-0.2 
-0.3 
-0.4 
-0.5 
-0.6 
-0.7 
-0.8 
-0.9 

x-location  [in] 

Figure  8  Total  deformation ,  time  step  =  5,  isometric  view  of  surface  plot 

Figure  9  is  similar  to  the  plots  above,  except  it  shows  the  cumulative  energy  required  to  perform  the  deformations  reported  in  the 
above  plots. 
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Figure  9  Total  energy  to  perform  deformation ,  time  step  =  5,  isometric  view 

Figure  1 1  shows  the  terrain  profile  surface  at  the  undisturbed  configuration  and  after  5  steps  are  performed.  Notice  that  the  vertical 
deformation  is  in  the  shape  of  an  ellipse  around  the  tire  contact  patch  when  using  the  subsoil  stresses  to  calculate  bulk  density  change 
and  sinkage,  rather  than  having  a  shape  identical  to  the  contact  patch,  as  is  the  case  if  the  tire  forces  are  applied  directly  to  the 
calculation  of  bulk  density  change  and  sinkage.  It  should  be  noted  that  the  propagation  of  subsoil  stresses  in  the  lateral  and  horizontal 
directions  is  reduced  when  a  larger  Frolich  parameters  is  used,  and  would  result  in  soil  deformations  that  are  more  condensed  around 
the  contact  patch. 

SUMMARY/CONCLUSIONS 

A  physics-based  tire  and  terrain  model  for  use  in  off-road  vehicle  dynamics  simulations  has  been  developed  and  presented.  The  tire 
model  is  a  simplified  mass-spring-damper  model,  while  the  soil  model  consists  of  multiple  components  that  model  the  vertical  and 
lateral  force/displacement  relationships.  The  soil  model  captures  the  change  in  soil  state  in  a  volumetric  sense,  and  does  not  make  the 
assumption  that  the  terrain  remains  homogenous  in  the  vertical  direction,  which  is  what  is  typically  done  in  most  empirical 
terramechanics  codes.  Repeated  loading  and  time-domain  effects  are  taken  into  account,  and  the  force/displacement  soil  equations  can 
be  utilized  to  calculate  the  required  energy  and  power  to  deform  the  terrain.  It  is  predicted  that  this  model  will  result  in  not  only  more 
realistic  mobility  simulations,  but  also  allow  for  a  novel  power  and  energy  analysis  technique  for  off-road  vehicles. 

Example  numerical  results  verify  that  the  soil  model  predicts  the  anticipated  visco-elasto-plastic  response  to  a  repeated  load,  and  that 
the  calculations  can  be  applied  to  a  terrain  database  to  produce  similar  results.  This  paper  represents  an  ongoing  effort  between  the 
University  of  Tennessee,  University  of  Wisconsin-Madison,  Mechsim,  and  the  US  Army  TARDEC  to  incorporate  the  physics-based 
VTIM  into  a  real-time  vehicle  dynamics  simulator.  As  this  is  a  work  in  progress,  updated  models  and  results  from  implementation  will 
be  reported  at  a  later  date.  Experiments  are  planned  to  measure  the  driveshaft  torque  of  a  vehicle  operating  on  deformable  terrain  as 
well  as  the  end  state  of  the  disturbed  terrain,  and  to  validate  the  proposed  models  by  comparing  the  experimental  and  simulation 
results  for  energy  and  power  required  to  perform  the  soil  deformation. 
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Figure  10  Cumulative  soil  deformation  [in],  steps  0-5,  (top  left  to  bottom  right) 


Figure  11  Undisturbed  and  final  terrain  profile 
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